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ABSTRACT 

We present a method to characterize statisticaUy the parameters of a detached binary sample - 
binary fraction, separation distribution, and mass ratio distribution - using noisy radial-velocity data 
with as few as two, randomly spaced, epochs per object. To do this, we analyze the distribution of 
ARVniax, the maximum radial- velocity difference between any two epochs for the same object. At low 
values, the core of this distribution is dominated by measurement errors, but for large enough samples 
there is a high-velocity tail that can effectively constrain the parameters of the binary population. 
Wc discuss our approach for the case of a population of detached white-dwarf (WD) binaries with 
separations that are decaying via gravitational wave emission. We derive analytic expressions for 
the present-day distribution of separations, integrated over the star-formation history of the Galaxy, 
for parametrized initial WD separation distributions at the end of the common-envelope phase. We 
use Monte Carlo techniques to produce grids of simulated ARVmax distributions with specific binary 
population parameters, and the same sampling cadences and radial velocity errors as the observations, 
and we compare them to the real ARVmax distribution to constrain the properties of the binary pop- 
ulation. We illustrate the sensitivity of the method to both the model and observational parameters. 
In the particular case of binary white dwarfs, every model population predicts a merger rate per star 
which can easily be compared to specific type-la supernova rates. In a companion paper, we apply 
the method to a sample of ^ 4000 WDs from the Sloan Digital Sky Survey. The binary fractions and 
separation distribution parameters allowed by the data indicate a rate of WD- WD mergers per unit 
stellar mass in the Galactic disk, ~ 1 x 10"^'^ mergers yr~^ Mg^, remarkably similar to the rate per 
unit mass of Type-la supernovae in Milky- Way-like galaxies. 

Subject headings: binaries: close, spectroscopic — white dwarfs — supernovae: general 



1. INTRODUCTION 

Stellar multiplicity is a fundamental piece in our cur- 
rent picture of stellar formation and evolution. Mod- 
ern studies of stellar multiplicity aim to constrain the 
fraction of stars with companions (or multiplicity frac- 
tion, /m), the distribution of separations, and the depen- 
dence of these parameters on variables like stellar mass, 
age, and metallicity. Different observational techniques 
are used to probe differ ent separation and flux contrast 
regimes (IDuauennov fc Mayor 1991; Malmrov_fc Kapl ani 
2005t iMason et all 120091 : iMetchev fc'HiUenbrandl f^Oa 
Raghavan et al.ll2010D . Short-period binary systems are 
of particular interest as the ancestors of future interact- 
ing binaries, from cataclysmic variables to novae, high 
and low-mass X-ray binaries, supersoft X-ray sources, 
and Type la supernovae. However, the fundamental 
properties of short-period binaries in the Galaxy are still 
poorly constrained. This has important implications 
for testing specific s cenarios for the formation of mul- 
tiple stellar systems (jTohlind 120021 : iGoodwin fc Kroupal 
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[2001 IBressert et all l2010l : IMarks et all l20Tll). stel- 
lar population synthesis models (IBruzual &: CharlotI 
[1991 ICOTlrov et all l2009t IMarks fc Kroupal 120 111) , and 
birth rate calculations for their interacting d escendants 
(jBelczvnski fc Taamll200l iRuiter et al.ll2009l ). 

The first modern spectr oscopic survey for close stellar 
binaries was performed by IDuauennov fc MavoH (fT99li ) , 
who examined 164 F and G-type stars, taking more 
than 4200 individual radial velocity (RV) measurements. 
Such surveys are extremely labor intensive, because they 
need to obtain enough RV measurements of each tar- 
get to either confidently discard multiplicity up to a 
certain period threshold, or to derive an orbital solu- 
tio n of sufficient quality. For example, the recent study 
bv lRaghavan et all ()2010[) examined 454 solar- type stars 
with different techni ques (including RVs), only a factor 
3 improvement over iDuquennov fc Mavod (|1991[ ) in al- 
most twenty years. The advent of large spectroscopic 
data bases like the Sloan Digital Sky Survey (SDSS, 
lYork et al.ll2000l) opens the possibility to take RV sur- 
veys to the next level, allowing millions of RV measure- 
ments of hundreds of thousands of different stars with 
well-calibrated and stable instrumental set-ups. 

In this paper, we present a method to characterize bi- 
nary populations statistically based on large samples of 
noisy RVs, with only a few epochs per target, but with- 
out followup observations, confirmations of real binaries, 
and derivations of orbital parameters for individual sys- 
tems. The observable that we analyze is the distribution 
of maximum RV differences, ARVmax, which is straight- 
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forward to obtain from a set of RV measurements for a 
sample of stars. This distribution contains information 
about the orbital velocity amplitudes of some members 
of some of the close binaries in a sample of stars. The use 
of velocity differences assures that the (uninteresting in 
the present context) systemic velocities of the stars are 
subtracted out of the dataset, along with any other time- 
constant velocity offsets (wavelength calibration errors, 
gravitational redshifts, etc.). Naturally, the more times a 
system is observed, the higher the chances of catching its 
full orbital RV variation amplitude, but even with only 
two epochs per system, some fraction of that amplitude 
will be probed. Apart from the observed RVs themselves 
and their sampling times, accurate knowledge of the dis- 
tribution of RV measurement errors is essential. We il- 
lustrate how the sensitivity of the method depends on 
t hese observational char a cterist ics. 

iTutukov fc Yungelsonl ()1986f ) outlined such an ap- 
proach schematically, as a way to discover close binary 
white dwarfs (WDs). The first application of this kind 
of technique (though i n a m ore rudimentary way) was 
by IMaxted fc Jeffried (|2005f ). A slightly different ap- 
pr oach that ha s been applied to SDSS data is described 
in iClark et al.l (|2011| ). Here, we present the method in 
the context of the problem of characterizing a population 
of detached WD binaries whose separations are decaying 
via gravitational wave emission. In a companion paper 
(Badenes & Maoz 2012; Paper II), we apply the method 
to a sample of ^ 4000 WDs from the SDSS with multiple- 
epoch spectra, we constrain the sample's binary popula- 
tion parameters, and we estimate its gravitational-wave- 
driven merger rate. 

2. MONTE CARLO SIMULATION OF THE 
ARVmaxDISTRIBUTION OF A BINARY POPULATION 

We now simulate the ARVmaxdistribution of an as- 
sumed binary population. Our methodology can be ap- 
plied to any binary population, but we will focus here on 
detached binary WDs in the Galactic disk, with the aim, 
in Paper II, of finding the regions of parameter space 
that describe a binary population that are allowed by 
the observed ARVmax distribution for a sample of SDSS 
WDs. We simulate WD binaries with properties drawn, 
in a Monte-Carlo process, from possible families of distri- 
butions of these properties, as detailed below. We then 
"observe" the simulated systems with the sampling se- 
quences and the velocity error distributions of the real 
data, to produce each model ARVmaxdistribution. A 
detached WD binary will merge, due to loss of energy 
and angular momentum to gravitational wave emission, 
within a time dictated by its separation and its compo- 
nent masses. A corollary of every population model will 
therefore be a prediction of the model's WD merger rate 
(whether sub-Chandrasekhar or super-Chandrasekhar), 
which can be directly compared to observed type-la su- 
pernova rates, or to the rates of other transient events 
that are candidates for the outcomes of such mergers. 
We describe below each step in this modeling process. 

Our modeling approach is distinct from that 
of "bin ary populatio n synt h esis" (BPS) calc u lation s 
such as [Tbonen et al.l (I2011D. iM enneke ns et all ()2010D . 
IWang et alj (|2010[) . or lRuiter et al. (200^. In BPS, one 
begins by simulating a population of main-sequence bina- 
ries, with a chosen mass and separation distribution, and 



one then attempts to follow the complex stages of stel- 
lar and binary evolution of each system, including mass 
transfer, mass loss, common envelope evolution, and so 
on. BPS calculations have many free parameters. Apart 
from the parameters that specify the initial conditions, 
there is a variety of parametrized ways to approximate 
the physics of various stages of evolution, particularly 
the common-envelope phases. Because of this variety 
and freedom, there is a great range among the predic- 
tions of different BPS calculations for the characteristics 
of the final WD populations. Instead, our approach is 
to parametrize in a simple mathematical way the prop- 
erties of the binary population at the end of its complex 
physical evolution - for binary WDs, at the end of the 
last common-envelope phase. Beyond that phase, there 
is only well-understood and easily modeled orbit decay 
due to gravitational wave emission. The general forms 
of our parametrizations for the component masses and 
separations at this evolved stage are guided by direct ob- 
servations, by BPS results, or by educated guesses. They 
are thus more realistic than those based solely on a par- 
ticular BPS realization, but they also allow investigating 
a larger parameter space for what the binary population 
might actually be like. The real RV measurements can 
then select the particular allowed regions of this param- 
eter space. 

2.1. WD primary mass 

For every simulated WD system (either single or bi- 
nary) we begin by assigning a primary mass ('primary' 
and 'secondary' refers here to the larger and smaller 
mass, respectively, not to which star will dominate the 
spectral energy distribution). We draw the primary 
mass, mi, fro m the observed distr ibution of WD masses 
determined bv lKepler et al.l (j2007[ ) for 1859 hot (effective 
temperature Toff > 12000 K) and bright {g < 19 mag) 
DA WDs in the DR4 SDSS catalog. We do not use the 
mass functions implied for the cooler WDs in each class, 
as Kepler et al. (2007) point out and discuss the uncer- 
tainty in the atmospheric modeling of those cool WD, 
which likely leads to a systematic over-estimate of their 
masses. The mass distribution is composed of four Gaus- 
sian components - a main, narrow, component centered 
at 0.58 M0 with la width of 0.047 Mq, and three weaker 
components centered at lower and higher masses. The 
latest version of the SDSS WD catalog, correspondi ng 
to D R7, has over 17000 entries (jKleinman et al.ll2009f n 
The iKepler et aTl (|2007D sample is a subset of the DR7 
WD sample that we actually analyze in Paper II, and so 
its mass distribution is likely quite representative of that 
of the WDs that we see in our sample. As a consistency 
check, we have measured the mass distribution for hot 
and bright DA WDs in the DR7 catalog, using the Toff 
and effective gravity, log 17, val ues fitted by Kleinman et 
al., and the cooling curves of iFontaine et al. (2001,) 0. 
We have confirmed with a KS test that this distr ibution 
is very close to the one obtained by iKepler et al.l for the 
DR4 WDs. 

* The catalog has not been published yet, but it was kindly 
made available to us by S. Kleinman (private communication). The 
version that we use here is from July 2010, and it has 17371 entries. 
^ We obtained an updated version o f these cooling curves from 
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It is hkely that, contrary to our assumption above, 
WD primaries do not have the same mass distribution 
as single WDs. Mass transfer and mass loss in close bi- 
naries can lead to either an increase or a decrease in the 
primary WD ' s fina l mass. In the BPS ca lculations of 
IClaevs et alj (|2011D . based on the code of llzzard et alj 
(|2006() . the primaries in WD pairs with separations of 
< 14i?Q have a broad mass distribution, between 0.3 
and 0.9 M0 (J. Claeys 2011, private communication). 
However, we find that even rather drastic changes in the 
primary mass function have only a small effect of the 
ARVinax distribution and on the merger rate. We have 
changed the center of the main Gaussian mass compo- 
nent to as low as 0.4 Mq or as high as 0.8 Mq , and find a 
negligible effect of the ARVmax distribution. Increasing 
the 1(7 width to 0.2 Mq also has only a small effect on 
ARVniax (at the level of, e.g., increasing the distribution 
by ~ 20% at ARVmax=300kms-^ see below). Our con- 
clusions will therefore not depend on the particular form 
we have assumed for the WD primary mass distribution. 
The weak dependence of ARVmax, in general, on the de- 
tails of the binary mass components, is discussed further 
in Section [3l 

2.2. Binary fraction 

One major binary population parameter to be con- 
strained by data is the fraction of objects of a class that 
is in binary systems. As we will see below, an analysis 
of a ARVniax distribution can be sensitive only to binary 
systems with velocities in the tail of the distribution, be- 
yond a "core" that is dominated by single systems and 
random velocity errors. For example, for the SDSS WD 
sample analyzed in Paper II, this threshold is at ARVmax 
>250 kms~^. For the range of possible binary com- 
ponent masses in a population, this velocity-difference 
threshold effectively puts an upper limit on the binary 
separations that can be probed. In an extreme-mass- 
ratio WD binary with masses mi = 1.2 Mq and 7712 = 0.2 
Mq, the secondary (7712) will achieve such orbital peak- 
to-trough velocity amplitude if the separation is ^0.05 
AU. For individual WDs with low-noise measurements 
and long temporal baselines, the SDSS data might be 
able to detect binaries with larger separations, but our 
analysis does not single out such systems. Therefore, we 
will define the binary population parameter, /bin, as the 
fraction of all WD systems (both single systems and bi- 
nary systems) that are binary systems with separations 
a < 0.05 AU. The length of an individual exposure deter- 
mines the minimum separation that we can detect. For 
the SDSS sample, the individual exposure times are ^15 
min. For the lowest WD masses, ~ 0.2 Mq, this cor- 
responds to separations of 10'^ km, or only about 10 
WD radii. Thus, the SDSS data can probe binary sep- 
arations ranging from close to Roche-lobe overflow, and 
out to 0.05 AU, which corresponds to periods of about 4 
days for typical WD masses. Since all double-degenerate 
(DD) systems that merge within a Hubble time have pe- 
riods shorter than ^ 12 hours, this is more than adequate 
to provide an estimate of the local WD merger rate. 

However, the binary fraction cannot be considered in- 
dependent of the masses of the WD components. From 
the initial-final mas s relation for stars and WDs (e.g. 
I Williams et aI1[2009l) . it is known that WDs with masses 
less than mum ~ 0.45 Mq have not had enough time 



to form in isolation over the age of the Universe, and 
therefore such WDs must be in binaries, with separa- 
tions that permit interactions in the course of the stellar 
evolution of the comp o nents . This was first shown di- 
rectly by iMarsh et al.l (jl995[ ). who found that five out 
of the seven WDs that they studied, with masses below 
0.45 Mq had WD comp a nions with periods P < 5 d. 
iRebassa-Mansergas et al.1 (|2011[ ) have shown that, in bi- 
naries composed of a main-sequence star and a WD of 
mass < 0.5 Mf:^, the period is ge nerally < a few days. 
Most recentlv. lBrown et al.l (|2012[ ) mined the SDSS pho- 
tometrically for WDs with masses < 0.25 Mq. After 
follow-up spectroscopy, they found that such low-mass 
WDs are always, or almost always, in close binaries, with 
periods of less than 1 day, and often with a relatively 
massive WD (or possibly neutron star) companion. The 
extre mely low-mass WDs actually found and followed 
up by iBrown et al.l (j2012| ) turn out to have even lower 
masses, always < 0.20 Mq. It is not known at what pre- 
cise mass the close (P < 1 d) binary fraction becomes 
100%. 

To account for this effect, we adopt a limit miim = 
0.25 Mq, such that when our simulated primary WD 
mass is smaller than m-nm, we always assign it to be in 
a binary. The fraction of the iKepler et al.l mass func- 
tion that is below this mass means that, effectively, we 
assign 0.07% of the WDs in the simulated sample to 
be in binaries in which one of the WDs is less massive 
than 0.25 Mq. Naturally, close companions likely exist 
also around all WDs with somewhat higher masses, e.g., 

< 0.35 Mq (which constitute 2.5% of the Kepler et al. 
2007 WDs) or < 0.45 Mq (9%). However, we do not 
know for a fact, for those higher masses, that the orbits 
are within the 0.05 AU limit that we have set above, when 
defining close binarity. Indeed, the high-ARVmax tail of 
our observed sample in Paper II is not dominated by low- 
mass objects, but rather by typical, ~ 0.6 Mq, WDs, as 
shown in FigHJ Thus, /bin in our simulations is the frac- 
tion of systems with separations a < 0.05 AU, and with 
both components with masses above mum ~ 0.25 Mq. 
Every simulation includes an additional population of 
extremely- low-mass (m < mum) WDs that are all in bi- 
naries and are not counted in /bin- To investigate how 
our results depend on the choice of miim, we have also 
calculated models with mum = 0.35 Mq and 0.45Mq. As 
was the case, above, with variations in the primary WD 
mass function, we find only weak effects on the ARVmax 
distribution, as long as the total fraction of binaries is 
equal (e.g., a model with mum = 0.25 Mq and /bin=0.05 
gives a very similar ARVmax distribution as a model with 
iTi\im = 0.35 Mq and /bin=0.025; for which the total 
binary fraction (including systems with a WD of mass 

< "T^iim) equals 0.05. Thus, our conclusions will not de- 
pend on the exact choice of mum. 

For the fraction /bin of the simulated WDs, as well 
as the extremely-low-mass WD binaries, we assign ad- 
ditional binary parameters, as described below. To the 
remaining 1— /bin WDs, we assign an orbital velocity of 
zero, and skip directly to the allocation of random ve- 
locity errors at several observing epochs (Section l2.6p . 
The maximum difference between these random errors 
for each such non-close-binary WD will then constitute 
the simulated ARVmax for that WD. 
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Fig. 1. — Mean and rms (shown with error bars) masses of the 
photometric primaries in the observed sample of Paper II, based on 
spectral modeling of the WDs, for various bins of ARVmax- The 
WDs with large ARVmax(> 200 km s~^) which reflect the close 
binaries in the sample, have masses ~ 0.6 Mq, similar to single 
WDs, whose mean and rms mass range is shown by the horizontal 
dashed and dotted lines. The bins without error bars are based 
on only one system each, and the 375 kms~^ bin is based on two 
systems. 

2.3. WD secondary mass 

The mass of the secondary, m2, is not hkely to be 
drawn from the same distribution as mi. Already in 
main-sequence binaries, it is clear that the secondary 
mass is not drawn from the initial mass frmction, but 
rather from a mass-ratio distrib ution that is approxi- 
mately flat (iRaghavan et al.ll201ClD . However, there is lit- 
tle in the way of observational or theoretical guidance for 
choosing the mass distribution of post-common-cnvclopc 
WD secondaries. Of the ^40 known DD systems, only 
10 are double-lined, i.e., with detected spectral features 
from both components. For these systems, both WD 
masses can be determined, but the number is still too 
small to say much a bout the mas s distr ibution. In the 
BPS calculations of iClaevs et al.l (|2011| ) the secondary 
WDs have a roughly flat mass-ratio distribution. We 
therefore choose the following parametrization. In cases 
where mi is above mii,„, if the simulated system is a bi- 
nary, we draw the secondary WD mass from a power-law 
distribution in mass ratio, 

P{q)(xq'^, q = 1712/1711, (1) 

with 1712 between mum and mi. The power-law index 
P is one of the parameters that we vary among the re- 
alizations of our simulation, in order to constrain the 
properties of the WD binary population. In cases where 
the primary in the Monte-Carlo draw was below mum 
(and therefore the star is always in a binary) , the second 
star is chosen with equal probability between 0.2 Mq 
and 1.2 Mq. In this scheme, since the typical mass pri- 
mary has a mass ~ 0.6 Mq, the secondary will have, 
on average, a mass of ~ 0.4 M©, reflecting the ob- 
serv ed com monnes s of low-mass W Ds in close binaries 
(e.g. iRebass a-M a.nsergas et al.l[2011| ) . 

For illustration purposes, we have plotted these pri- 
mary and secondary mass distributions in Figure ^ 
with P ~ 0. We have overlaid the standard bound- 
aries between He, CO, and ONe cores in isolated WDs 
(0.5 and 1.0 Mq, respectively), but we note that these 
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Fig. 2. — Probability density distribution of primary (blue) and 
secondary (red) masses in our model binary WD population. A 
binary fraction of /ijin=0.05 is assumed in this example, and this 
sets the height of the extremely low-mass WD peak (< 0.25 M©), 
in which WDs are always in close binaries. A flat mass-ratio dis- 
tribution ((9 = 0) is also assumed in this example. We note that 
the two distributions are not independent. For illustration, the 
theoretical boundaries between He, CO, and ONe WDs in isolated 
stars are shown with vertical lines. 

boundarie s might no t apply to close binary systems (see 
[Prada Moroni fc St ranicro, 2009.) . 

2.4. Separation distribution 

Next, we assign to each simulated DD system a sep- 
aration, and hence we need to consider what are the 
possible separation distributions for close binary WDs. 
The separation distribution at the time the WDs emerge 
from their final common envelope phase is unknown ob- 
servationally, while theoretically it is a longstand ing, 
complex, and unsolved problem (see IIvanovaTl20 1 ll for 
a recent review). The close WDs may undergo eithe r 
one or two common-envelope phases ([Woods et al.ll20l3) , 
which could, in principle, lead to a complicated separa- 
tion distribution. Nevertheless, BPS cal culations, as well 
as so me more sophisticated treatments ([Delove fc TaamI 
|2010[ ). suggest a power-law separation distribution with 
a negative index, over the range of separations that we 
consider here , ~ 0.1 — 10 Rq. In the BPS calculations of 
IClaevs et al.l (l20Tll) . over the range of separations that 
we consider, the post-common-envelope initial WD sep- 
arations indeed are roughly constant per logarithmic in- 
terval (J. Claeys, private communica t ion) . The same is 
true in the BPS models of lYungelsonI (j2010t) , at least for 
separations above ^ 1 R© (L. Yungelson, private com- 
munication). The ^ t~^ Type-la supernova delay-time 
distributions generally predicted by BPS models for the 
DD channel would not arise if the WD initial separation 
distribut ions were not approxima tely of the above form 
(see, e.g. lMaoz &: Mannuccill20lH) . 

Therefore, we will assume an initial WD separation 
distribution that is a power law, with an index that is a 
free parameter to be constrained by observations. What- 
ever the initial distribution, orbital decay due to gravi- 
tational wave emission will immediately begin to modify 
it, as all of the orbits shrink and the innermost systems 
merge. Furthermore, the actual distribution at any par- 
ticular time will be the sum of the distributions of many 
populations of different ages, that have evolved over dif- 
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ferent amounts of time. We now calculate the form of 
this eveolved, time-integrated, distribution. 

2.4.1. Evolution of a binary separation distribution due to 
gravitational wave losses 

The separation a of two point masses, mi and TO2i 
in a circular-orbit binary, will shrink over time due to 
gravitational wave energy loss as 

da _ K 256 



g, — r 5 'tiiTO2(toi + m-2), (2) 



where G is the gravitational cons tant and c is the speed 
of light (iPeters fc Mathewill963D . From integration, the 
time t for the system to evolve from separation a' to 
separation a obeys 

a'4 -a^ = Kt. (3) 

Suppose a population of WD binaries is formed at a time 
t = (this could be, e.g., a population of WD binaries 
in the Galaxy that have just emerged from the common 
envelope phase). The population has an initial distribu- 
tion of separations n'(a'). For simplicity, we will assume 
the initial distribution of WD masses is independent of 
separation. Systems with separations in the range a' to 
a' + da' will migrate, after a time <, to a bin a to a -I- da 
in the evolved distribution n(a,t). Conservation of the 
number of systems (except for those systems that reach 
a = and merge) requires that 



n(a,t)da ~ n'{a')da' , 



or 



n(a,t)=n'(a')^=n'(a')(^)' 
= n'[{a^ 



da 



(4) 



(5) 



(6) 



(a4 + Kt)^^ ■ 
If the initial distribution is a power law, 

n'(a') (X a'", 
then 

n(a,t) DC a3(a4 + if t)("-3)/4_ (7) 

The evolved distribution at time t is thus approximately 
a broken power law. For separations a ^ {Kty^^ (i.e., 
much larger than those that can merge within time t) , the 
distribution will have approximately the original power- 
law slope, n{a) ~ a". At separations a <^ (Kt)^/'^, on the 
other hand, n(a) ~ a^ (see left panel of Figure [3]). The 
merger rate versus time from this single-age population 
will be controlled by the short separation pairs for which 



dn dn da 

— = — — oc n[a,t) a 
dt da dt 



-3 ^ t("-3)/4^ 



(8) 



This is a well-known result for the delay-time distri- 
bution of the gravitational-wave-induced mergers of a 
single-age population, having an initial separ ation dis- 
tribution that is a power law of index a (e.g. iGreggiol 
[200llTotam et al. 1120081) . 

Let us consider now a series of binary WD populations, 
each with an initial separation distribution n'(a') oc a'", 
being produced at a rate R{t) between t — and the 



present age of the Galaxy, to- The present-day distribu- 
tion will be 



N{a) = / R{to - t)n{a,t)dt 
Jo 



oc / R{to-t)a^{a^ + Kt)^"-^^^^dt. 



(9) 



Assuming, for instance, a constant star-formation rate 
over the age of the Galaxy, then also R{t) — const. (Even 
if the star-formation history is "bumpy" , the WD forma- 
tion history will be the convolution of the star-formation 
history with a broad, ^ kernel, that describes the 

WD supply [e.g., Pritchet et al. 2008], and which will 
smooth out the WD production rate). The integral then 
gives 

7V(a;)ocx4+«[(l + a;-4)("+i)/4-l], a ^ -1, (10) 



or 



where 



N{x) oc ln(l -t- x-"), 



(11) 



(12) 



is the separation in units of the separation of binaries 
that will merge within the age of the Galaxy. For exam- 
ple, for to = 10 Gyr and mi = m2 = 0.55 Mq, x = 1 
corresponds to oq = 0.01 AU, or 1.5 x 10^ km, or about 
150 WD radii. The right panel on Figure [3] shows N{x) 
for various values of a. N{x) is, again, approximately a 
broken power law, with index a ai x I. Atx^l the 
power- law index is 3 for a > —1, and a + 4 for a < — 1. 

2.4.2. Choice of binary separation 



We use the functional forms in Eqns. [TOlITT] to model 
the possible present-day distributions of DD separations, 
for various indices a of the initial power-law distribu- 
tions at formation. In realizations of our simulation, we 
draw the separation of specific component binary masses 
from the present-day distribution for a particular value 
of a, with a between a^m = 2 x 10^ km (contact) and 
flmax = 0.05 AU. For the purpose of producing simu- 
lated radial velocities, binaries with periods < 15 min 
are assigned zero orbital velocity, as the SDSS exposure 
length prevents detection of velocity differences in such 
cases. A practical consideration in this process is the 
large dynamic range that the distribution P{a) can as- 
sume over this range in a. As a result, very few simu- 
lated binaries might be assigned separations in ranges 
that have low probability, and the model expectation 
values for the the velocity differences due to those bi- 
naries will have large Poisson errors. To avoid this, we 
populate the distribution evenly with simulated systems 
among four decades of separation a (i.e., from Omin to 
flmax/lOOO, from flniax/lOOO to aniax/100, ctc). Each bi- 
nary system is given a relative weight according to the 
integral of the separation distribution over the decade it 
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Fig. 3. — Left: Example of the gravitational-wave-driven time evolution of the separation distribution of a population with an initial 
power-law separation distribution of index a = —2.0 (Eq. O. In this example, both WDs are of 0.55 Mq. Right: Separa tion distributions 
for various values of a, time-integrated over a constant star-formation rate over 10 Gyr in the Galactic disk (Eqns. llOllll l. Separations are 
relative to the initial separation of two 0.55 MqWDs that will merge within 10 Gyr, a = 0.01 AU. 



Given two masses and an orbital separation, Kepler's 
law gives the period 



T = 27r 



G{mi + ?7i2) 



1/2 



and the circular orbital velocities, 



■"1 



27ra TO2 

T TOl -I- 7712 ' 



V2 



2'Ka mi 
T mi + m2 



(13) 



(14) 



We assume circular velocities for simplicity, but also 
because this is the expectation for close binaries that 
have likely undergone circularization by tidal forces and 
common-envelope evolution. For example, there are es- 
sentialy no binaries with P < 12 d that h ave any appre- 
ciable eccentricity ('Raghavan et al.' ' 20101 ) . It has been 
recently suggested (Thompson 2010j) that there might 
be a preponderance of triple systems amon g binary WDs 
that would induce elliptical orbits via thelKozaJ ([l96l 
mechanism, leading to a large population of systems with 
short merger times. We defer the exploration of this sce- 
nario to future work. Equation [3] with a = gives the 
merger time, tmoigoj of the simulated system. The merger 
rate per WD in the simulated sample is obtained by not- 
ing, for a given set of parameters, the fraction of all of the 
systems in the simulation that merge within a set time 
window, divided by that time. Each system is weighted 
according to its decade in separation (see section [2.4.21 
above). A separate tally is conducted to calculate the 
merger rate of only those binaries whose summed masses 
exceed the Chandrasekhar mass, which may be of rele- 
vance for the type-la supernova rate. For a > — 1, the 
merger rate is approximately constant, and therefore the 
time window for numerically calculating the rate is arbi- 
trary, as long as it is shorter than to- The constancy of 
the rate can be seen from Eq. [TOl by noting that mergers 
come from systems with x < 1, for which N{x) ^ x^, 
and the merger rate is 



dN dN da , o 

— = — — cx N{a) a"'^ 
dt da dt ^ ' 



const. (15) 
For a < — 1, the merger rate falls with time, but quite 



slowly for values of a that are not too steep. 



— N{a) a~ 

dt ^ ' 



t(a+l)/4 



(16) 



so an accurate numerical merger rate is still obtained in 
the above scheme. 

The merger rate for a given combination of /bin and a 
can also be roughly estimated analytically. If a > — 1, 
the majority of pairs have large separations, with long 
times until merger. The merger rate is therefore set by 
the integrated effect of old systems. All binaries with 
separations of a < ao at their formation time will con- 
tribute (at a constant rate, as we saw) to the current 
merger rate, where oq is the maximum separation bi- 
nary that merges within tg- Every choice of component 
masses gives a different value of oq. However, for the 
range of possible component masses, ag is between 0.005 
and 0.018 AU and a value of oo = 0.01 AU is typical. 
The typical time until merger of each system is ig- The 
merger rate per observed WD is therefore 



1 dN ^, /bin 

TVwd dt (1 — /bin)^o la'^^ n{a)da 



(17) 



/b 



(ao/flmin) 



(1 — fhin)to (flmax/amin) 



Q-l-1 



1' 



for n(a) = a". If /bin is small, and (as in the present 
case), Cniax and ag are both much larger than amin, then 
the merger rate is roughly 



1 dN 
N^~dt 



/bill 



ag 



a+1 



(18) 



This gives values in good agreement with the merger 
rates from the numerical Monte Carlo calculation. One 
can also see from Eq. [T8|that, in a plot of the parameter 
space of a vs. log /bin, curves of constant merger rate 
will be straight lines with slope of 1/ log(ag/aniax) ~ 1-4, 
for ag = 0.01 AU and a,nax = 0.05 AU (see Paper II). For 
a ^ — 1 ; most WD binaries are formed with small sepa- 
rations and therefore merge within a time much shorter 
than to- In this case, the merger rate will be controlled 
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by the supply rate of new WDs. which in turn is set by 
the star-formation rate. 

2.6. Inclination, photometric primary, temporal 
sampling, and velocity error 

We next apply observational effects to each simulated 
binary sytem. First, a line-of-sight inclination i of the 
perpendicular to the orbital plane is chosen from a dis- 
tribution 

P{i) cx sini, (19) 

and the line-of-sight velocity is reduced by sini. 

We then need to decide which of the two WDs will be 
the photometric primary. This could be either the less 
massive WD, because it has larger surface area and/or it 
is younger and hence hotter; or the more massive WD, 
because it cools more slowly due to its small surface area 
and is hence hotter. In practice, these effects compete 
against each other, and it is difficult to determine which 
of the two components will dominate the spectral en- 
ergy distribution. From an observational point of view, 
the distribution of absolute magnitudes as a function of 
WD mass in the DR7 SDSS catalog shows a large spread 
about the mean at all masses, although low mass WDs 
(below ~0.35 Mq) do seem to be intrinsically brighter 
(see Figure IH) . From a theoretical p oint of view, the cool- 
ing curves of iFontaine et al.l ()2001f ) also predict that, in 
coeval DD systems, WDs with masses below ^ 0.35 Mq 
will remain substantially brighter than their more mas- 
sive counterparts for several hundred Myr (in the SDSS 
g filter - the effect is enhanced in r, and diminished in 
u). After about 1 Gyr, the slower cooling of more mas- 
sive WDs takes over and makes them brighter, but by 
this time the predicted magnitudes become fainter than 
the cutoff in our samples for most of the volume probed 
by SDSS. To summarize, it seems reasonable to assume 
that low-mass WDs, if present, will have a tendency to 
dominate the spectral energy distribution of DD systems, 
but in other cases either of the components may be dom- 
inant. In our Monte Carlo runs, we therefore make the 
less massive WD the photometric primary when its mass 
is below 0.35 Mq, but decide randomly between the two 
WDs when the less massive WD is above this limit. 

Once the photometric primary has been selected, we 
sample its line-of-sight velocity with the actual distribu- 
tion of temporal samplings in the SDSS data. We do 
this by choosing at random a particular observation pat- 
tern (number of epochs and time between epochs) from 
the sample, with a random phase assigned to the first 
epoch of the sinusoidal RV curve. Figure [S] shows the 
distributions of the number of epochs and of the tempo- 
ral baselines per object for the SDSS sample of Paper II. 
To each simulated velocity measurement, we add a ran- 
dom error that we draw from a Gaussian distribution, 
with the variance of the Gaussian drawn from the ac- 
tual distribution of measurement errors for the observed 
sample (see figure 1 in Paper II for this distribution for 
the SDSS sample). Finally, for every simulated system 
or non-binary WD, we find the minimum and maximum 
observed velocities and calculate ARVmax- 

For every parameter combination that defines a binary 
WD population model, we produce 10^ WD systems (ei- 
ther single or binary, according to /bin), and find the frac- 
tional prediction for each bin in the model ARVmax dis- 




0.6 0.8 
WD Mass [M 

su.. 

Fig. 4. — Distribution of absolute g magnitudes as a function 
of mass for the hot DA WDs in the DR7 catalog. We only con- 
sider, in this figure, objects with T^ff above 12000 K and g mag 
brighter than 19, for which spectroscopic masses can be modeled 
reliably. The overlaid histogram shows the mode (most frequent 
g magnitude observed) and the standard deviation in bins of 0.1 
Mq. There is a large spread in absolute magnitudes at all masses, 
but low mass WDs (below ~ 0.35 Mq) do seem to be intrinsically 
brighter. 

tribution. Multiplied by the observed WD sample size, 
this gives the expectation value for that bin. 

3. RESULTS 

Figure [6] shows the simulated ARVmax distribution for 
a model binary population with a particular set of pa- 
rameters (/bin, a, /?), when a sample of that population 
is sampled with a particular set of empirical temporal se- 
quences and a particular velocity error distribution. The 
observational parameters chosen are those of the SDSS 
sample presented in Paper II, and the model parame- 
ters are one set among those that reproduce the data. 
Along with this distribution, we plot the predictions for 
a model with no binaries in it. We see that the mod- 
eled ARVmax distribution consists of two parts. At low 
ARVmax, there is a "core" that is dominated by random 
velocity errors, that have been applied to systems that 
are single, or that have low ARVmax because of one or 
more causes (low orbital velocity, low inclination, inop- 
portune sampling). Beyond this core, the distribution 
has a "tail" that reveals those close binaries that hap- 
pened to have large orbital velocities, favorable inclina- 
tions, and temporal sampling that caught their velocity 
variations. For any observational sample, one can always 
calculate this zero-binaries core. When compared to the 
observed ARVmax distribution, it immediately reveals if 
and where in the data there exists a tail of real binary 
systems. A statistical comparison between the data and 
a grid of model distributions can then select the regions 
of binary population parameter space that are allowed or 
ruled out by the data. 

A reliable estimate of the RV errors is essential, as oth- 
erwise underestimated errors can masquerade as a tail, or 
overestimated errors can hide a real binary population. 
This is also illustrated in Fig. IH where we show several 
no-binary predictions that use incorrect error distribu- 
tions, i.e., error distributions that are different from the 
one used to generate the distribution with binaries. 

The temporal sampling density of the survey will nat- 
urally affect both the core and the tail of the ARVmax 
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Fig. 5. — Temporal sampling characteristics for the DR7 WD sample analyzed in Paper II, and used in the example simulations here. 
Left: Distribution on number of epochs per object. Most WDs have only two or three epochs, but there is a tail of objects with more 
epochs. Right: Distribution of maximum time differences between epochs. A 45 min interval is most common, but longer timescales, of 
order an hour or of several days, are probed in many cases. 
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Fig. 6. — Simulated distributions of ARV max. Solid line shows 
the distribution resulting for a WD population with binary fraction 
/bin=0.05, separation distribution power-law index a = —1, and 
mass-ratio distribution power-law index /3 = 0, when this distri- 
bution is observed with the temporal sequences and RV errors of 
the DR,7 sample presented in Paper II. The dashed line shows the 
"core" of the distribution, obtained by setting /bin=0 (no binaries, 
except for the extremely low-mass WDs, which are always in bina- 
ries), and reflects the part of the full distribution that is due to RV 
errors alone. Colored curves show the core distributions obtained 
when assuming incorrect error distributions (various narrow distri- 
butions of la error ranges, as marked). Accurate characterization 
of the errors is thus essential for characterizing the binary popula- 
tion by means of the signal in the tail of the ARVmax distribution. 

distribution. The more epochs per object, the greater 
the chance of catching the fuh RV variation range of 
each system. The core will also grow, but only as the 
square root of the number of epochs. Since the fully re- 
vealed RV range reaches saturation beyond some number 
of epochs, while the core ARVmax continues to grow as 
random errors are added in quadrature, there will be a 
limiting typical number of epochs, N, beyond which the 
technique is no longer efficient for characterizing the pop- 
ulation, and one can, instead, attempt to fit RV curves to 
each object. This will happen when \/N<7^ ^ ARVmax -i-i 
where the latter is the highest value of ARVmax for which 
a sample has some exemplars, and (t„ is the typical ve- 
locity error. Fig [7] illustrates how the core and the tail 
of the distribution change when, instead of the full sam- 
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Fig. 7. — Dependence of ARV max distribution on temporal sam- 
pling. Blue curves show the full ARVmax distribution (solid curves) 
and the /bin=0 core (dashed curve), for a chosen set of binary pop- 
ulation parameters, and with the actual DR7 Paper II sample with 
its temporal sampling characteristics (same as Fig. |6]) Red curves 
show the distribution and core when every object is observed on 
two epochs only. 

pling of the SDSS dataset, every object is sampled only 
twice. 

Figure [8] shows how the ARVmaxdistribution depends 
on the binary population parameters /bin, ct, and /?. 
Qualitatively, increasing /bin and decreasing a both 
have the same effect of increasing the number of small- 
separation binaries, and therefore of raising the high 
ARVmax tail of the distribution. Quantitatively, if we 
were dealing with a population of binaries with single 
values of component masses, single values of inclination 
angle, numerous sampling per object, and no measure- 
ment errors, then the ARVmax distribution tail would 
behave as 

dN dN da ^^^^ 



dv 



da dv 



and the broken power-law separation distribution, N{a), 
would lead to a broken power-law orbital velocity distri- 
bution. Recalling that, for a > — 1, the small-separation 
part of the separation distribution behaves as N(a) oc a^, 
and taking the Keplcrian v cx a~^/^, we would then ex- 
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Fig. 8. — Dependence of the ARV max distribution on binary pop- 
ulation parameters. Top: Dependence of binary fraction, /bin. The 
blaclc solid curve is for no binaries, except for those that include a 
< 0.25 MqWD. The dashed curve is with no binaries at all, and 
shows the bare "core" of the distribution, due to the radial veloc- 
ity errors. Middle: Dependence on separation distribution power 
law index, a. Bottom: The weak dependence on the mass-ratio 
distribution's power-law index, /3. 

pact the tail to fall steeply as ^ (for all a > —1), 
with its amplitude depending on /bin/(a + 1) (the latter 
factor entering through the normalization of the distribu- 
tion). The a" branch of N{a) at a > oq would transform 
to a t;~2Q-3 power law at velocities v < vq, where vq is 
the orbital velocity of a binary with separation oq. In re- 
ality, the distribution of component masses (which leads 



to a range of values for ap), the projection of velocities 
due to inclination, and the sampling (both of which ef- 
fects move binaries in the distribution to lower values of 
ARVinax), and the velocity errors (which again mix the 
distribution), will lead to a ARVmax tail with a slope 
that behaves differently than above, and does depend on 
a. For a < -1, N{a) cx a°'+^, and the ARVmax dis- 
tribution tail falls more gently, as ^ ■y-2(Q-fii)^ ^^le 
idealized case. Thus, in principle, a can be discerned 
in data with low-enough velocity errors, and with large 
enough samples, such that the slope of the tail can be 
measured accurately. 

As seen in Fig. [U the ARV,„ax distribution is weakly 
dependent on /3, the power-law index of the binary mass- 
ratio distribution. For a given choice of primary mass, 
the secondary's velocity will depend on the mass ratio q 
as (1 + q)"^/^. For the typical mi < 0.75 M0 primary 
mass, the secondary mass is constrained to the range 
from mi down to mum = 0.25 M0, so q is between about 
1/3 and 1. Going from high positive values of /3 that 
favor q near 1, to very negative /3's that concentrate q 
for all binaries to be near 1/3, will induce an increase of 
■\/3/2 = 1.22 in the secondary's velocity V2, and an even 
smaller relative decrease, by 0.93, in vi, if it is the pri- 
mary that is selected as the photometric primary. Thus 
when changing between the extreme values of (3, about 
half of the binaries will have their ARVmax increase by 
22%. and half will decrease by 7%, or a net increase by 
7% in the ARVmax of each bin in the distribution. Even 
for the highest velocities that we consider, this is a small 
change, comparable to the typical velocity errors. Fur- 
thermore, as we saw above, the ARVmax distribution is 
roughly a power law. The transformation v ^ v' = kv, 
where fc is a constant, does not affect the shape of power- 
law distribution in v. 

4. CONCLUSIONS 

We have used Monte Carlo simulations to study how 
the distribution of maximal radial velocity differences, 
ARVmax, can characterize a close binary population, in 
a radial velocity survey where a large number of stars 
are observed, but with only a few epochs per star, and 
with potentially large RV errors. Our focus has been on 
a survey for double WDs of the kind that we analyze 
in Paper II, using the SDSS DR7 WD sample, which 
has served here as our example survey in terms of obser- 
vational parameters. As part of this modeling process, 
we required a realistic representation of the present day 
separation distribution of close WD binaries whose sepa- 
rations have evolved due to gravitational wave emission. 
We have therefore derived analytical expressions for the 
separation distribution of a population of WD binaries 
that is continuously formed, its orbits decay, and some of 
its members merge, over the lifetime of the Galaxy. With 
these simulations, we have determined how the ARVmax 
distribution depends on the binary population, which 
we have characterized using three parameters (describing 
close-binary fraction, separation distribution, and mass- 
ratio distribution), and on the observational parameters 
of the sample - RV errors and temporal sampling pat- 
tern. 

Our main findings are: 

1. The ARVmax distribution has a core region, that 
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is produced by the random RV errors, and a tail, 
that can reveal the presence of some of the close 
binaries in a sample. The power of the technique 
is that, even with large errors and only few (or 
even just two) epochs per object, the close binary 
population can reveal itself in the tail, provided 
that the number of objects in the sample is large 
enough. This tail permits statistically constraining 
the properties of the population, even with sparse 
and noisy data, and without detailed followup work 
on candidate binary systems. 

2. Accurate knowledge of the distribution of radial 
velocity errors is essential to model properly the 
ARVniax distribution, and thus deduce the contri- 
bution of the real binary population to the tail. 

3. Steep distributions in initial separation (very neg- 
ative a, with most binaries at small separations) 
and populations with a large close-binary fraction 
(large /bin) will both result in an increase in the 
amplitude of the ARVmax tail. There may thus be 
some degeneracy in the determination of these two 
parameters. A change in a, however, also changes 
the shape of the tail, and hence, in high-quality 
datasets (many objects, small errors), these param- 
eters may still be individually constrained. 

4. The ARVmax distribution depends weakly on the 
details of the distributions of the component binary 
WD masses - the distributions of primary masses 
and of secondary mass ratios, which in any case 
are chosen from a a relatively small dynamic range. 
This binary population characteristic thus cannot 
be constrained by this kind of survey data. Con- 
versely, not knowing the distribution of mass ratios 
does not affect adversely our ability to constrain 
the other binary population parameters. 

5. The ARVmax approach allows to estimate the 
merger rate of a close binary population, based on 
noisy RV data with few epochs, of the kind we con- 
sider. This is possible without follow-up observa- 
tions to obtain full binary parameter solutions for 
candidates, and without necessarily even finding a 
single binary that will merge within a Hubble time. 
This ability is a result of the statistical nature of 
our approach. 

In Paper II we apply the methods presented here to 
the observed SDSS DR7 WD sample, wc set constraints 
on the local population parameters of close- WD binaries 
in the Galaxy, and we derive the merger rate of this pop- 
ulation, both in general and for particular component 
and total mass ranges. We show that the local rate of 
WD mergers with super-Chandrasekhar total masses is 
an order of magnitude lower than the local Type-la su- 
pernova rate. The local merger rate of all WDs, however, 
is remarkably similar to the Type-la suparnova rate. A 
large fraction of all WD mergers are between CO and CO 
WDs. with total masses not far from the Chandrasekhar 
limit. If sub- Chandrasekhar mergers result in a Type-la 
supernova explosion, we may have identified their domi- 
nant progenitors. 



Apart from Type-la supernova explosions, other pos- 
sible outcomes of WD mergers can be tested with 
our methodology - R Corona Borealis stars (e.g. 
iLongland et al.ir2011| ). or highly magnetic WDs, which 
have been postulated to result from WD mergers (e.g. 
iGarcfa-Berro et"aI1l2012l and references therein). Some 
7 ± 3% of local WDs have magnetic fields greater than 
lO'' G (jKawka et al.ll2007|) . Assuming these magnetic 
fields decay on very long timescales, WD mergers pro- 
ducing all of this population would need to occur, over 
10 Gyr, the age of the Galaxy, at a rate of ~ (7 ± 3) x 
10^^^ yr^^ per WD. In Paper II, our fit to the observed 
ARVmax distribution for the WDs in SDSS indicates a 
WD merger rate, for total merged masses of < 1 M0, 
of l^Q y X 10"^^ yr^i per WD. The rate may thus be 
sufficient to explain at least some, and perhaps even all, 
magnetic WDs with such mergers. 

As another example, about half of WD s with masses 
below 0.45 Mp, appear t o be si ngle (iMaxted et al.l 
2000t 'Napiwotzki et all [20071: iKilic et al.l I2007D . 
Nelemans & Tauris (19981) have suggested formation of 
such low-mass single WDs via interaction between a 
giant star and a close massive-planet or brown-dwarf 
companion, follow ed by evaporation or tidal disruption 
of the companion. IKilic et all ()2007[) have proposed that 
these WDs have evolved from metal-rich stars whose 
evolution was truncated by severe mass los s on the giant 
branch. Alternatively, llben et all ()1997[ ) have raised 
the possibility that the single low-mass WDs are the 
merger products of even-lower WDs. Our constraints 
on merger r a tes c an test this last scenario. In the 
iKepler et all ()2007D WD mass function, about 8% of 
the WDs are in a Gaussian component that peaks at 
around 0.4 M0. As suming tha t half of these WDs are 
single, then in the llben et all ()1997f ) scenario, about 
4% of the WD population would be the outcome of 
very-low-mass WD mergers. This is a similar fraction 
to the one invoked above in the case of magnetic WDs, 
and therefore would require a similar merger rate. Our 
calculations, however, show that the rate of mergers 
with total masses in the range 0.3-0.5 Mq is four orders 
of magnitude below the one required by this scenario. 
This comes about because each of the merging WDs 
would need to be below 0.3 Mq, and such WDs are rare. 
The longer gravitational orbit decay time at these low 
masses further lowers the rate. One could circumvent 
this argument by in voking large mass loss during the 
merger process (e.g. iFrver et al.l l2010( ). so that more 
massive and com mon WDs could be involved. However, 
iDan et all ()2012l) find negligible mass loss in their 3D 
hydrodynamic simulations of WD mergers. 

The tools that we have introduced here can also easily 
be adapted to the characterization of stellar multiplicity 
based on large RV surveys in other settings. Exam ples 
are ongoin g surveys like LAM OST (|Su et al.lll998l) and 
APOGEE fPrieto et aIll2008D. a nd planned ones, such 
as BigBOSS (Schlegel et al.ll20ll . 
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